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Abstract 
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complex  cases  with  multiple  layers  can  complicate  the  numerical  implementation.  A  control  volume 
analysis  is  used  in  the  present  work  to  derive  a  finite  volume  approach  to  the  1-D  thermal  analysis 
problem  using  generalized  coordinates.  Details  of  the  numerical  implementation  are  also  discussed  to 
highlight  how  this  technique  can  simplify  the  overall  solution  process.  Finally,  two  example  problems 
are  used  to  validate  the  accuracy  of  the  approach,  one  of  which  involves  analysis  of  a  multilayer  TPS 
problem  with  non-uniform  grid  spacing. 
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Nomenclature 


a  =  thermal  diffusivity 

e  =  emissivity 

p  =  density 

a  =  Stefan-Boltzman  constant 

6  =  time  integration  weighting  factor 

£  =  spatial  coordinate,  computational  plane 

cp  =  heat  capacity 

E  =  energy  rate  with  respect  to  time 

p  =  directional  heat  flux  vector 

Fo  =  Fourier  number 

k  =  thermal  conductivity 

n  =  surface  normal  vector 

q  =  heat  flux  (per  unit  area) 

qg  =  heat  generation  rate  (per  unit  volume) 

S  =  surface  area 

T  =  temperature 

t  =  time 

u  =  generic  property  in  a  control  volume 

V  =  volume 

x  =  spatial  coordinate,  physical  plane 


Subscripts/Superscripts 

cond  =  conduction 

conv  =  convection 

g  =  internal  generation 

j  =  grid  index 

L  =  left  side 

n  =  time  index 

R  =  right  side 

rad  =  radiation 

x  =  derivative  with  respect  to  spatial  coordinate  x 
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1.  Introduction 


The  solution  of  the  transient  one-dimensional  heat  diffusion  equation  is  an  elementary 
problem  taught  in  many  courses  on  heat  transfer.  While  exact  solutions  are  possible  for  a  subset  of 
problems,  engineering  applications  typically  involve  using  numerical  techniques  to  obtain  an 
approximate  solution  to  the  heat  equation.  For  a  simple  analysis  involving  a  single  material  with 
uniform  grid  spacing,  finite  difference  methods  are  straightforward.  For  more  complex  engineering 
problems  where  the  computational  grid  spacing  is  not  uniform,  and/or  where  different  materials  are 
being  considered  (i.e.  composite  slabs),  it  will  be  shown  that  the  implementation  of  a  finite  volume 
technique  using  a  generalized  (i.e.  curvilinear)  coordinate  transformation  can  help  simplify  the 
numerical  approach. 

Although  several  authors  have  published  material  discussing  the  numerical  solution  of  the 
heat  diffusion  equation  using  finite  volume  techniques  [1]  or  two-dimensional  generalized 
coordinates  [2],  and  industry  standard  thermal  analysis  tools  such  as  FIAT  [3]  use  a  finite  volume 
solution  technique,  the  author  is  unaware  of  a  reference  describing  the  derivation  and  finite  volume 
numerical  implementation  of  the  transient  one-dimensional  heat  diffusion  equation  using  generalized 
coordinates.  The  present  work  is  therefore  aimed  at  addressing  this,  as  well  as  presenting  an  example 
of  an  engineering  application  of  the  technique  for  multi-material  composite  slabs. 


2.  Finite  Volume  Formulation  Using  Generalized  Coordinates 

A  finite  volume  approach  is  an  intuitive  technique  that  can  be  associated  with  the  control 
volume  analysis  that  can  be  used  to  derive  the  heat  diffusion  equation.  In  this  approach,  the  quantities 
of  interest  in  a  given  control  volume  are  typically  defined  as  a  de  facto  averaged  value.  For  example, 
consider  a  variable  u(x,  y,  z)  that  varies  as  a  function  of  location  within  a  three-dimensional  control 
volume.  The  “averaged”  value  of  u  can  be  found  by  integrating  the  function  over  the  entire  volume: 


udV  =  uV  (1) 

Jv 

where  iZ  represents  the  averaged  value  of  u. 

To  apply  this  approach  to  the  heat  diffusion  equation,  first  consider  an  energy  rate  balance  for 
a  control  volume,  which  can  be  written  assuming  conservation  of  energy  as 


(2) 


In  equation  (2),  the  rate  of  energy  stored  in  the  control  volume  (. Estored )  is  equal  to  the  net  energy 
rates  into  the  control  volume  due  to  conduction  (Econd),  convection  (Econv),  radation  (Erad),  and 
internal  energy  generation  (Eg)-  The  rate  of  energy  stored  can  be  written  as 


d[  [  dT  ( dT\ 

Estored  =  Yt  Jy  PCP  TdV  =  Jy  pCP~^dV  =  [pCP  ~^)V 


(3) 
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where  V  represents  the  size  of  the  control  volume.  The  term  pcp(dT / dt)  is  the  change  with  respect 
to  time  of  thermal  energy  in  the  material  per  unit  volume,  with  the  overbars  on  the  variables  denoting 
“averaged”  values  for  the  control  volume  as  defined  previously.  Note  that  in  equation  (3),  it  is 
assumed  that  the  material  density  and  specific  heat  do  not  change  over  the  time  interval  being 
evaluated. 

The  net  energy  rates  due  to  conduction,  convection,  and  radiation  can  be  represented  by 
integrating  the  net  heat  flux  entering  the  control  volume  due  to  the  appropriate  physical  phenomena. 
For  example,  the  net  energy  rate  due  to  conduction  can  be  written  as 

Econd  ~  I  Fcond  '  n  dS  (4) 

Js 

where  Fcond  represents  the  directional  heat  flux  due  to  conduction,  n  is  the  outward  pointing  surface 
normal,  and  S  represents  the  total  surface  area  of  the  control  volume.  Similar  equations  can  be  written 
for  the  convective  and  radiative  net  energy  rates. 

Finally,  the  net  change  in  energy  due  to  internal  heat  generation  can  be  written  as  the 
averaged  value  of  the  heat  generation  rate  per  unit  volume  multiplied  by  the  control  volume  size: 


Jv 


qg  dV  =  qgV 


(5) 


Using  equations  (3)-(5),  the  conservation  law  defined  in  equation  (2)  can  be  rewritten  in 
integral  form  for  the  control  volume  as 


j  PCp  TdV  j  Fcond  tr  dS  +  j  Fconv  ■  n  dS  +  j  Frad  ■  n  dS  +  j  qg  dV 


(6) 


In  the  case  of  a  one-dimensional  control  volume,  the  volume  reduces  to  a  length  Ax,  and  the 
surface  area  reduces  to  unity;  therefore  equation  (6)  can  then  be  rewritten  as 


-  dT\  r 
Pcp  ^)Ax  = 


condin  F condout )  +  ( FConvin  FConvout 


)  +  ( Fradin  Fradout)  +  (7) 


where  the  subscripts  in  and  out  refer  to  the  conductive,  convective,  and  radiative  heat  fluxes  into  and 
out  of  the  control  volume  in  the  x-direction. 

Using  Fourier’s  law  for  conductive  heat  transfer,  equation  (7)  can  then  be  rewritten  as 


+  qg  Ax 

where  the  derivative  dT / dx  in  equation  (8)  is  evaluated  at  the  control  volume  interfaces. 


(8) 
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Non-uniform  spacing  of  the  finite  volume  cells  can  complicate  the  implementation  of 
numerical  approaches  to  solving  Equation  (8).  Therefore,  a  coordinate  transformation  can  be  applied 
using  a  generalized  coordinate  (x)  that  is  a  known  function  of  x;  typically  the  transformation  is 
defined  such  that  <f  (x)  has  uniform  spacing  to  simplify  the  numerical  analysis.  Applying  the  chain 
rule  of  differentiation,  Equation  (8)  can  be  converted  to  generalized  coordinates  as 


7  d dT\  /  d%  dT\ 
V  k  dx  df)  in  V  ^  dx  3<f  / 

+  {fradin  ~  F radout )  +  Rg^x 


out  j 


+  ( FconVin  FConvouO 


(9) 


3.  Numerical  Implementation 

The  previous  section  discussed  the  derivation  of  the  one-dimensional  heat  diffusion  equation 
using  a  control  volume  approach.  From  a  numerical  standpoint,  this  can  be  thought  of  as  equivalent  to 
a  cell-centered  finite  volume  approach,  where  the  center  of  each  cell  is  assigned  an  integer  index 
value  j,  and  the  cell  interfaces  have  a  fractional  (e.g.  j  +  1/2)  index. 


In  the  solution  of  equation  (9),  the  coordinate  transformation  is  defined  such  that  f  (x)  has 
uniform  spacing  of  A<f  =  1.  The  grid  transformation  metrics  in  equation  (9)  can  be  numerically 
approximated  at  the  cell  center  as 


Ey.  =  — 
Sx)  dx 


Cy+i/2  tj- 1/2  _  1  _  1 

-*7  +  1/2  —  Xj- 1/2  -*7  +  1/2  —  Xj- 1/2  AXy 


(10) 


Likewise,  the  grid  metrics  at  the  cell  interfaces  can  approximated  numerically  in  a  similar  fashion: 

fj-fj- 1  1 


_  d% 

^t-l/2  -  fa 


j~  1/2 


Xj  —  Xj Xy  —  Xy_! 


_  d% 

j+1 /2  dx 


fy+i  -  0 


(ID 


/+1/2  X/+1 


Xy  Xy  +  1  -  Xy 


For  the  remainder  of  the  discussion  in  this  section,  the  convective  and  radiative  heat  fluxes 
are  assumed  to  be  zero,  as  these  are  typically  incorporated  into  the  boundary  conditions  for  many 
engineering  problems  of  interest.  Additionally,  the  internal  heat  generation  term  is  assumed  to  be  zero 
(i.e.  no  heat  sources  exist  in  the  material).  Equation  (9)  can  then  be  rewritten  for  a  computational  cell 
with  index  j  as 


(  dT\ 

17  dT\ 

/  dT\ 

AXy  = 

j 

-1 

7-1/2 

1 

|rc> 

,  * 

1 

7  +  1/2- 

(12) 


The  time  derivative  can  be  approximated  numerically  using  a  first  order  backwards 
difference;  the  conductive  heat  fluxes  at  the  cell  interfaces  can  be  constructed  numerically  based  on 
the  temperature  in  the  cells  on  each  side  of  the  interface.  This  allows  equation  (12)  to  be  written  as 
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(13) 


Pjcp 


rptl  _rptl-l 

j 


At 


A  Xj 


=  (1-0)  k  if.  i 


'-rTl—'L _ 'T'TL—  1 

D'+i  V 


7+9 


■/+2 


+  9  (  kj+\^j+\ 


rp7 l  'T’TL 

‘j+ 1  D 

Af 


"k.  l^X; 


li  li- 1 

A^ 


'T'Tl _ 'T’TL 

li  1 

A/ 


where  n  refers  to  the  current  time  step  in  a  numerical  integration,  and  n  —  1  corresponds  to  the 
previous  time  step.  The  weighting  factor  0  is  defined  between  0  and  1,  where  0  =  0  represents  an 
explicit  time  integration  approach,  and  0  =  1  is  fully  implicit. 

Note  that  the  material  density  and  heat  capacity  in  equation  (13)  are  defined  at  the  cell  center 
j,  whereas  the  thermal  conductivity  is  defined  at  the  cell  interfaces  j  ±  1/2.  One  may  be  tempted  to 
use  an  arithmetic  mean  to  determine  the  thermal  conductivity  at  the  cell  interfaces.  Patankar  [1], 
however,  has  shown  that  this  can  lead  to  incorrect  results  when  analyzing  composite  slabs  with  abrupt 
changes  in  the  thermal  conductivity  between  materials,  and  that  a  better  approach  is  to  consider  a 
control  volume  between  the  cell  centers  of  the  two  cells  adjoining  the  interface  as  a  composite  slab, 
and  apply  a  steady-state  one-dimensional  heat  conduction  analysis.  Assuming  that  the  cell  interfaces 
are  midway  between  the  cell  centers,  this  analysis  simplifies  such  that  the  interface  thermal 
conductivity  is  the  harmonic  mean  of  the  thermal  conductivities  in  the  two  cells  adjoining  the 
interface,  which  can  be  written  as 


Xj±1/2  ~ 


(14) 


Due  to  the  fact  that  the  coordinate  transformation  is  defined  such  that  A<f  =  1,  and  noting 
from  equation  (10)  that  Ax;-  =  1/ /  for  a  given  cell  with  index  j,  rearranging  and  simplifying  terms 
allows  equation  (13)  to  be  rewritten  for  an  explicit  integration  scheme  with  0  =  0  as 


™ n  _  rptl—  1 

I  4  ; 


+  ■ 


At 


Cxy  kj+  l/2f; 


f ’pn—1  _T,n—l 
Xj+l/2Vj+l  ‘j 


)  -  /2(7rI-7Ti1  j 


(15) 


The  explicit  integration  scheme  shown  in  equation  (14)  is  a  simple  function  of  the 
temperature  distribution  at  the  previous  time  step,  but  the  Fourier  number  must  be  less  than  or  equal 
to  0.5  (for  a  one-dimensional  analysis)  to  maintain  numerical  stability  [4],  where  the  Fourier  number 
is  defined  as 


a  (At) 
(Ax)2 


¥At)  (r  \2 

PiCPj  ^ 


(16) 


Depending  on  the  specific  problem  being  analyzed,  the  thermal  diffusivity  may  be  a  function 
of  temperature  and  location  in  a  multi-material  composite  slab,  and  the  spacing  in  the  x-direction  may 
not  be  uniform.  Therefore,  for  an  explicit  analysis,  the  stability  constraint  should  be  evaluated  on  a 
per-cell,  per-time  step  basis  to  determine  the  maximum  time  step  possible  in  the  integration  of 
equation  (15). 
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For  an  implicit  time  integration  technique  with  0  =  1,  equation  (13)  can  similarly  be 
simplified  to 


Note  that  equation  (17)  is  a  linear  system  of  equations  that  must  be  solved  simultaneously  to  update 
the  temperature  for  time  step  n,  but  has  no  numerical  stability  limitations  on  the  time  step  size. 
Equation  (17)  results  in  a  tridiagonal  matrix  of  coefficients  for  the  system  of  equations,  which  can  be 
solved  efficiently  using  standard  inversion  techniques  such  as  the  Thomas  algorithm. 


4.  Calculation  of  Cell  Interface  Temperatures  and  Boundary  Conditions 

While  the  finite  volume  approach  discussed  above  makes  intuitive  use  of  control  volume 
analysis  to  derive  a  numerical  solution  technique,  one  challenging  aspect  lies  in  the  way  the  boundary 
conditions  are  implemented.  As  the  material  being  analyzed  is  discretized  into  finite  volume  cells  to 
solve  equation  (15)  or  (17),  the  cell-centered  approach  discussed  means  that  the  temperature  is 
updated  at  the  cell  centers.  The  question  then  arises  as  to  how  to  determine  the  temperature  at  the 
interfaces  between  cells,  as  the  edges  of  the  material  (which  are  cell  interfaces)  are  typically  prime 
areas  of  interest  in  an  engineering  analysis. 

Similar  to  the  discussion  of  calculating  the  cell  interface  thermal  conductivity,  the  use  of 
steady-state  one-dimensional  heat  conduction  is  useful.  Consider  the  control  volume  shown  in  Figure 
1,  where  the  boundaries  of  the  control  volume  correspond  to  locations  halfway  between  a  finite 
volume  cell  interface  and  the  centers  of  the  two  finite  volume  cells  adjoining  the  interface.  A  steady- 
state  one-dimensional  heat  conduction  equation  can  be  written  for  this  control  volume  similar  to 
equation  (12),  except  that  the  left  hand  side  of  the  equation  is  set  to  zero  due  to  the  steady-state 
assumption,  and  the  fluxes  are  constructed  numerically  taking  into  account  that  the  control  volume  is 
half  the  size  of  a  normal  finite  volume  cell: 

(  dT\  (  dT\ 

\k^xdVL~\k^xdVR~°  (18) 
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One-Dimensional  Control  Volume 


Cells  7+V 2 


Interface 

Figure  1.  Depiction  of  a  one-dimensional  control  volume  used  to  determine  the  temperature  at 
the  cell  interfaces. 


Again  approximating  the  derivatives  numerically,  this  reduces  to 


Tj+i  Tj+ 1/2 


Tj+ 1/2  Tj 


As  the  L  and  R  boundaries  of  the  control  volume  lie  within  finite  volume  cells,  the  thermal 
conductivity  at  the  control  volume  boundaries  is  assumed  to  be  the  “averaged”  value  in  each  cell  as 
defined  in  equation  (1).  Therefore,  equation  (19)  can  be  rewritten  (dividing  out  the  constant  Af/2  and 
approximating  the  grid  metrics  at  the  L  and  R  boundaries  numerically)  as 

kj+i%xj+1{Tj+i  ~  Tj+ 1/2)  ~  kjl;Xj(Tj+l/2  -  Tj)  =  0  (20) 

Finally,  rearranging  terms  then  allows  the  cell  interface  temperature  to  be  calculated  as 

kj  +  l^X  T  +  l^jf  +  l  T  kjfX  :Tj 

W=^f - +k(  '  (2D 

For  TPS-related  engineering  applications,  the  problem  is  typically  such  that  a  multilayer  TPS 
material  with  boundaries  on  each  end  is  subject  to  convective  heat  fluxes  (e.g.  due  to  a  high 
temperature  flow  over  the  material)  and  radiative  heat  fluxes  (e.g.  when  the  surface  temperature  is 
significant).  To  apply  these  types  of  boundary  conditions  to  the  solution  of  equation  (15)  or  (17),  a 
straightforward  technique  is  to  use  “ghost”  cells  that  extend  past  the  boundary,  as  shown  in  Figure  2. 
For  this  example,  a  steady-state  heat  flux  balance  can  be  written  at  the  boundary  wall  (index  j  = 
1/2) as 


RconcL  —  Rconv  T  IJracL 
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Ghost  Cell 


o 
•  II 

• 

7=1 

<N 
•  II 

/ 

^  / 

\ 

7  =  1/2  7=3/2 

Wall  Boundary  Cell  Interface 

Figure  2.  Example  of  a  one-dimensional  ghost  cell  used  to  apply  a  wall  boundary  condition 


Using  Fourier’s  law  and  the  Stefan-Boltzmann  law,  this  can  be  rewritten  as 

/  dT\ 

[~k^x  gf)  =  Qconv  ~  (23) 

where  the  variables  in  equation  (23)  are  evaluated  at  the  boundary  wall.  Note  that  while  the 
convective  heat  flux  is  typically  written  in  terms  of  the  heat  transfer  coefficient  and  difference 
between  the  recovery  and  wall  temperatures,  for  the  purposes  of  the  present  work  it  is  treated  as  a 
“general”  heat  flux  input  into  the  surface  that  is  not  dependent  on  wall  temperature.  Likewise,  for 
simplicity  the  temperature  of  the  surroundings  used  to  determine  the  radiative  heat  flux  is  also 
assumed  to  be  zero  in  equation  (23).  Note  that  while  the  steady-state  heat  flux  balance  used  above 
only  considers  boundary  condition  heat  fluxes  due  to  conduction,  convection,  and  radiation,  other 
thermodynamic  phenomena  (e.g.  ablation  heating  or  cooling)  could  be  incorporated  into  equation  (22) 
as  appropriate. 


Using  the  indices  listed  in  Figure  2,  equation  (23)  can  be  rewritten  (approximating  the  spatial 
derivative  numerically  and  using  subscripts  to  denote  the  grid  indices)  as 


'ki/2^x1/2(^i  W  —  qc 


’ 1/2 


g£1/2(T\/2) 


(24) 


Rearranging  terms  allows  the  temperature  in  the  ghost  cell  to  be  calculated  as 


T0  =  T,+- 


Qc 


’’1/2 


—  as 


1/2 


(W 


kl/2%X 


(25) 


1/2 


The  ghost  cell  at  j  =  0  is  essentially  an  “imaginary”  mirror  image  cell  used  to  allow  the  use 
of  equations  (15)  or  (17)  with  no  special  treatment  at  the  boundaries  (aside  from  normal  handling  of 
boundary  condition  values).  Therefore  it  is  assumed  to  be  the  same  size  and  have  the  same  material 
properties  as  the  cell  at  j  =  1,  i.e.  the  thermal  conductivity,  emissivity,  and  grid  metrics  at  the  wall 
are  identical  to  the  averaged  properties  for  the  cell  at  j  =  1  (e.g.  kj=1/2  =  /c;=0  =  kj=1).  Equation 
(21)  can  then  be  used  to  calculate  the  temperature  at  the  wall  as 


(26) 


_  +  k&Jo  _T1  +  T0 

1/2  k^Xl  +  k^Xl  2 

In  equation  (25)  the  temperature  at  T0  is  a  function  of  the  temperature  at  7l/2  due  to  the 
radiative  heat  flux,  yet  in  equation  (26)  the  temperature  at  T1/2  is  a  function  of  T0.  Therefore,  in  the 
numerical  application  of  the  boundary  condition  where  the  radiative  heat  flux  at  the  wall  is  non-zero, 
equations  (25)  and  (26)  need  to  be  solved  using  an  iterative  technique.  While  in  some  cases  equations 
(25)  and  (26)  can  be  successively  iterated  until  the  values  of  T0  and  Tl/2  reach  a  constant  value  within 
a  specified  tolerance,  convergence  may  be  slow  or  difficult  to  achieve  in  situations  with  large 
temperature  gradients  and  high  surface  temperatures.  The  use  of  a  robust  root  finding  technique  (such 
as  Brent’s  algorithm  [5])  can  eliminate  this  issue. 


The  same  technique  discussed  above  can  be  used  to  derive  the  boundary  condition  at  the 
j  =  jmax  +  1/2  boundary  wall  of  the  grid  using  a  ghost  cell  at  j  =  jmax  +  1.  This  results  in 


Tjmax+l  ~  Tjmax  T 


Qconvjmax+ 1/2  a£jmaxVjmax+ 1 


:(j'jmax+ 1/2) 


k 


■jmax 


+l/2$x 


(27) 


jmax+1/2 


and 


rp  _  Tjmax+1  T  Tjmax  rjQ) 

f  jmax+1/2  —  2 

which  again  must  be  solved  using  an  iterative  technique  for  boundaries  with  non-zero  radiative  heat 
flux  until  T jmax  and  Tjmax+1/2  reach  a  constant  value. 

Note  that  for  an  implicit  time  integration  technique,  equations  (25)  and  (27)  supplement  the 
linear  system  described  in  equation  (17)  to  provide  equations  for  the  temperature  at  the  boundaries. 


5.  Validation  and  Example  Cases 


To  assess  the  validity  of  the  solution  technique  described  above,  a  thermal  analysis  tool  was 
created  in  Excel  using  Visual  Basic  for  Applications.  A  simple  validation  case  was  run  considering  a 
problem  in  which  a  semi-infinite  slab  of  copper  (, k  =  401  W/m*K,  a  =1.17x10  4  m2/s)  at  an  initial 
uniform  temperature  of  Tinitiai  =20  °C  is  exposed  to  a  constant  net  heat  flux  input  of  qnet  =  3xl05 
W/m2  on  one  side.  This  problem  is  identical  to  Example  5.8  described  in  Ref.  [4],  which  notes  that 
the  analytical  solution  for  T  (x,  f )  at  a  given  time  t  and  location  x  from  the  wall  can  be  written  as 


T(pc,  t) 


'  initial 


+ 


7  *7 11  e  t  \/  n  t  /  77" 


-exp 


Qnetx 

k 


erfckh) 


(29) 


The  finite  volume  method  discussed  in  the  present  work  was  used  to  solve  the  one¬ 
dimensional  unsteady  heat  diffusion  equation  in  the  copper  slab,  where  qnet  is  assumed  to  be 
equivalent  to  qCOnv  +  Rrad  in  equation  (22).  For  the  purposes  of  this  analysis,  the  results  from  the 
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explicit  algorithm  shown  in  equation  (15)  were  compared  to  the  analytical  solution  defined  above,  as 
well  as  to  the  explicit  finite-difference  algorithm  results  listed  in  Ref.  [4],  This  comparison  was 
performed  for  the  case  where  Fo= 0.25  and  Ax=75  mm,  which  results  in  At=12  s.  The  results  obtained 
by  solving  equation  (15)  for  the  first  9  cell  interfaces  of  the  finite  volume  cells  nearest  the  wall  are 
identical  to  that  listed  in  Ref.  [4]  for  the  finite  difference  solution  on  9  nodes  (located  at  the  same  x 
locations  as  the  finite  volume  cell  interfaces).  Finally,  Figure  3  plots  the  temperature  distribution  in 
the  material  at  t=120  s  as  a  function  of  distance  from  the  wall.  As  can  be  observed,  the  finite  volume 
solution  is  equal  to  the  finite  difference  solution,  and  both  are  reasonable  approximations  to  the 
analytical  results  obtained  by  equation  (29). 


Distance  from  wall,  m 


Figure  3.  Temperature  distribution  in  copper  slab  example  at  time  t=l  20  s. 


While  the  validation  case  above  demonstrates  that  the  finite  volume  formulation  can  produce 
results  equivalent  to  those  generated  using  a  node-centered  finite  difference  approach,  the  primary 
benefit  of  the  technique  described  in  the  present  work  is  in  the  application  to  more  complex 
engineering  problems  involving  multi-material  composite  slabs  and  non-uniform  grid  spacing. 
Consider  the  example  of  a  gliding  reentry  vehicle  with  a  thermal  protection  system  (TPS) 
configuration  consisting  of  the  materials  shown  in  Table  1.  A  thermal  analysis  can  be  conducted  to 
assess  the  temperatures  in  the  material  as  the  reentry  vehicle  is  subjected  to  the  heat  flux  indicated  in 
Figure  4.  For  this  purpose,  the  implicit  version  of  the  finite  volume  technique  described  in  the  present 
work  was  used;  the  same  example  was  also  run  with  the  industry  standard  thermal  analysis  program 
FIAT  [3]  to  confirm  the  results.  Both  analyses  used  the  same  material  properties,  outer  boundary  heat 
flux  input  show  in  Figure  4,  and  an  adiabatic  inner  boundary  condition.  At  the  outer  boundary,  the 
radiative  heat  flux  was  also  included  (which  required  an  iterative  solution  of  (25)  and  (26)  as 
discussed  in  Section  4). 
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Table  1.  Materials  in  TPS  Composite  Slab  Example 


Layer 

Material 

Thickness,  cm 

Outer  Layer 

Carbon-Carbon 

1.27 

Interior  Layer 

Carbon  Felt 

0.95 

Interior  Layer 

Min-K  Insulation 

4.57 

Inner  Layer 

Titanium 

1.27 

Time,s 

Figure  4.  Heat  flux  input  to  outer  boundary  for  TPS  composite  slab  example. 


The  outer  and  inner  boundary  temperature  time  history  for  this  TPS  composite  slab  example 
is  compared  in  Figure  5;  it  can  be  seen  that  the  results  from  FIAT  are  essentially  identical  to  those  of 
the  present  work.  (For  the  FIAT  analysis,  the  default  option  of  having  FIAT  generate  the 
computational  grid  was  used;  this  resulted  in  64  grid  points  for  FIAT  vs.  the  23  cells  used  in  the  finite 
volume  results  shown  in  Figure  5.)  The  temperature  profile  through  the  composite  slab  is  also  shown 
at  several  different  times  in  Figure  6;  the  points  on  each  “Finite  Volume”  line  correspond  to  the  cell 
interfaces  used  in  the  calculation  (23  cells  in  total),  and  demonstrate  how  the  technique  can  produce 
accurate  results  while  allowing  the  use  of  non-constant  cell  spacing  with  relative  ease.  (Note  that  for 
simplicity  the  results  for  this  TPS  composite  slab  example  assume  a  non-ablating/non-catalytic 
surface,  which  may  not  be  realistic  for  the  materials  and  temperatures  in  this  example.  Flence  these 
results  should  be  taken  as  a  first-order  thermal  assessment,  with  additional  modeling  required  to 
assess  the  effects  due  to  thermochemical  processes.) 
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Figure  5.  Thermal  analysis  results  for  TPS  composite  slab  example. 
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Figure  6.  Temperature  distribution  within  TPS  composite  slab  example  at  different  times. 
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6.  Conclusions 


The  present  work  has  outlined  a  finite  volume  approach  to  solving  the  transient  one¬ 
dimensional  heat  diffusion  equation  using  generalized  coordinates.  As  demonstrated  above,  it  can  be 
seen  how  this  technique  can  simplify  the  numerical  solution  of  multi-material  TPS-related 
engineering  problems  using  non-uniform  grid  spacing.  This  can  be  useful  in  allowing  the  number  of 
grid  cells  and  cell  spacing  in  each  layer  of  a  multi-material  TPS  analysis  to  be  tailored  to  the  accuracy 
required,  with  no  need  to  maintain  constant  spacing  throughout  the  material. 

While  beyond  the  scope  of  the  present  work,  it  should  be  noted  that  one  of  the  advantages  of 
deriving  the  finite  volume  solution  approach  from  a  control  volume  analysis  is  that  the  extension  to 
two-  and  three-dimensional  problems,  as  well  as  problems  with  non-hexahedral  cells,  is  relatively 
intuitive.  As  equation  (6)  implies,  the  problem  reduces  to  considering  the  rate  of  energy  stored  in  the 
cell  and  the  heat  transfer  across  each  of  the  cell  interfaces.  The  present  work  considered  heat  fluxes 
due  to  conduction,  convection,  and  radiation,  but  other  heat  fluxes  (e.g.,  due  to  ablation,  surface 
catalycity)  could  be  added  in  a  similar  fashion.  While  additional  dimensions  introduce  added 
complexity  in  the  grid  metrics,  such  as  the  multi-dimensional  application  of  the  chain  rule  of 
differentiation  and  the  need  to  model  fluxes  in  the  three  local  coordinate  directions  for  each  cell  face 
for  a  three-dimensional  cell,  standard  techniques  used  to  solve  multi-dimensional  fluid  mechanics 
problems  on  generalized  coordinates  are  applicable  [6]. 
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